
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import LineString, Point,Polygon,MultiLineString
from shapely.ops import nearest_points, substring,linemerge
import math
import rtree
from shapely import hausdorff_distance
from polylabel import polylabel
from shapely import minimum_bounding_radius,segmentize
from scipy.spatial import ConvexHull

""" 初匹配 """
################################################################################################################
#############################################################

def gen_matchrelation(gdf1, gdf2, buffer_size,arr,overlap_threshold=0.9):
    """找到gdf2中满足条件的line1，即gdf2：line1a，line1b,line1c...."""
    ##构建缓冲区
    buffer1=custom_buffer(gdf1, buffer_size)
    buffer2=custom_buffer(gdf2, buffer_size)

    # 构建树
    buffer1_idx = build_rtree_index(buffer1)
    buffer2 = buffer2.sort_index()
    match_list = []
    #遍历寻找
    for row in buffer2.itertuples():
        overlaps = buffer1_idx.intersection(row.geometry.bounds)
        lines_list = []
        for idx in overlaps:
            feature1 = buffer1.loc[idx]
            if row.geometry.intersects(feature1.geometry):
                line1 = gdf1.loc[gdf1['myid'] == feature1.myid].iloc[0]
                line2 = gdf2.loc[gdf2['myid'] == row.myid].iloc[0]
                overlap_lenth1 = line1.geometry.intersection(row.geometry).length
                overlap_lenth2 = line2.geometry.intersection(feature1.geometry).length
                # 匹配条件
                if (overlap_lenth1 / line1.geometry.length) > overlap_threshold or (overlap_lenth2 / line2.geometry.length) > overlap_threshold:
                    if  attribute_comparison(line1, line2,arr)==1:
                        lines_list.append(feature1.myid)
        if lines_list:
            lines_list.insert(0, row.myid)
        if lines_list:
            match_list.append(lines_list)

    df_lst = [pd.DataFrame([row]) for row in match_list]
    match_df = pd.concat(df_lst, axis=0, ignore_index=True)
    return match_df

""" 构建空间查询树 """
def build_rtree_index(gdf):
    bounds = [geom.bounds for geom in gdf.geometry]
    idx = rtree.index.Index((i, b, None) for i, b in enumerate(bounds))
    return idx

""" 属性对比 """
def attribute_comparison(line1, line2,arr):
    a=0
    if line1[arr]==line2[arr]:
        a=1
    return(a)
def calculate_max_buffer(geom):
    """最大缓冲半径计算"""
    if geom.area > 20:
        points = np.array(geom.exterior.coords).astype(np.float64)
        (_, radius) = polylabel([points], with_distance=True)
        return radius * 0.8
    else:
        return (geom.area / (minimum_bounding_radius(geom) * 4)) * 0.8

def custom_buffer(gdf, buffer_size):
    """改进的缓冲区生成方法"""
    # 复制原始数据框
    buffered_gdf = gdf.copy()
    
    # 定义缓冲逻辑
    def _buffer_logic(geom):
        if geom.geom_type == 'LineString' and geom.is_closed:
            try:
                polygon = Polygon(geom.coords)
                if polygon.is_valid:
                    max_radius = calculate_max_buffer(polygon)
                    if buffer_size > max_radius:
                        return polygon.buffer(buffer_size, cap_style=3)
            except:
                pass
        return geom.buffer(buffer_size, cap_style=2)
    
    # 应用缓冲并修复拓扑
    buffered_gdf.geometry = buffered_gdf.geometry.apply(_buffer_logic)
    buffered_gdf.geometry = buffered_gdf.buffer(0)  # 修复拓扑
    
    return buffered_gdf

def process_df(match_df):
    """ 处理匹配关系的df """
    # 将NaN值替换为空字符串
    match_df.fillna('', inplace=True)  
    # 转置，去除重复行，再次转置
    match_df = match_df.T.drop_duplicates().T

    new_match_df = pd.DataFrame(columns=["col1", "col2"])
    # 遍历每一行
    for i, row in match_df.iterrows():
        col1 = row[0]
        for j in range(1, len(row)):
            # 获取第j列数据
            colj = row[j]
            if colj:  # 如果colj不为空
                new_row = {"col1": col1, "col2": colj}
                # 将新行添加到new_match_df中
                new_match_df = pd.concat([new_match_df, pd.DataFrame([new_row])], ignore_index=True)

    return new_match_df

##################################################################################################################
###############################################################################################################################################################################

""" 开放式计算方法 """

def cal_semiangle(point1,point2):
    # 输入起点和终点计算范围(0-180)
    angle = math.atan2(point1[1]-point2[1], point1[0]-point2[0])
    angle=math.degrees(angle)
    if angle<=0:
        angle += 180
    return angle

def ge_shapecode(line):
    """输入为geopandas的GeoSeries对象"""
    line_geo = line.geometry
    MRrect = line_geo.minimum_rotated_rectangle
    
    # 计算外包矩形尺寸和圆心
    if isinstance(MRrect, LineString):
        x_diff = MRrect.coords[0][0] - MRrect.coords[-1][0]
        y_diff = MRrect.coords[0][1] - MRrect.coords[-1][1]
        or_angle = int(cal_semiangle( MRrect.coords[0],MRrect.coords[-1]))
    elif isinstance(MRrect, Point):
        # print("错误ID",line.myid)
        x_diff, y_diff = 1, 1  # 异常处理
        or_angle = 90
    else:
        corners = MRrect.exterior.coords
        x_diff = Point(corners[1]).distance(Point(corners[2]))
        y_diff = Point(corners[1]).distance(Point(corners[0]))
        if x_diff>y_diff:
            or_angle = int(cal_semiangle(corners[1],corners[2]))
        else:
            or_angle = int(cal_semiangle(corners[1],corners[0]))

    R=math.sqrt(x_diff*x_diff+y_diff*y_diff)/2
    o = Point(MRrect.centroid.x, MRrect.centroid.y)

    # 绘制同心环区域
    radii = [0.2 * R * i for i in range(1, 6)]  # 10层同心环
    width = 0.2 * R
    angle_inter =60

    # radii = [0.1 * R * i for i in range(1, 11)]  # 10层同心环
    # width = 0.1 * R
    # angle_inter =20

    # radii = [0.0625 * R * i for i in range(1, 17)]  
    # width = 0.0625 * R
    # angle_inter =12

    # radii = [0.05 * R * i for i in range(1, 21)]  
    # width = 0.05 * R
    # angle_inter =6


    # 初始化向量和扇形容器
    vector = [0] * int((360/angle_inter))
    sectors = []

    for level, radius in enumerate(radii, start=1):
        inner_radius = radius - width
        for angle in range(or_angle, or_angle + 360, angle_inter):
            # 生成外弧的边界点
            outer_arc_points = [
                (o.x + radius * np.cos(np.radians(angle)), o.y + radius * np.sin(np.radians(angle))),
                (o.x + radius * np.cos(np.radians(angle + angle_inter)), o.y + radius * np.sin(np.radians(angle + angle_inter)))
            ]
            # 生成内弧的边界点（逆序）
            inner_arc_points = [
                (o.x + inner_radius * np.cos(np.radians(angle + angle_inter)), o.y + inner_radius * np.sin(np.radians(angle + angle_inter))),
                (o.x + inner_radius * np.cos(np.radians(angle)), o.y + inner_radius * np.sin(np.radians(angle)))
            ]
            ring_sector = Polygon(outer_arc_points + inner_arc_points)
            # sectors.append(ring_sector)
            if ring_sector.intersects(line_geo):
                # 计算角度区间索引（0-11）
                sector_idx = ((angle-or_angle) // angle_inter)
                # 仅保留最大层级（若当前层级更高）
                if vector[sector_idx] < level:
                   vector[sector_idx] = level
    return vector, sectors, MRrect, o  

def cosine_sim(line_a, line_b):
    vec1, _, _, _=ge_shapecode(line_a)
    vec2, _, _, _=ge_shapecode(line_b)
    vec1 = np.array(vec1)
    vec2 = np.array(vec2)

    norm_vec1 = np.linalg.norm(vec1)
    norm_vec2 = np.linalg.norm(vec2)
    
    # 计算余弦相似度
    if norm_vec1 == 0 or norm_vec2 == 0:
        return 0  # 避免除以零的情况
    cosine_sim = np.dot(vec1, vec2) / (norm_vec1 * norm_vec2)
    # 长度比值
    if line_a.geometry.length != 0 and line_b.geometry.length != 0:
        cosine_sim=cosine_sim*min(line_a.geometry.length/line_b.geometry.length,line_b.geometry.length/line_a.geometry.length)
    else:
        cosine_sim=0

    return cosine_sim

################################################################################################################################################################################

""" 开放(截断)-闭式匹配。参数：容差大小和缓冲区距离 """


def merging_line(contour_gdf, line_b0, line_a0):
    """ # 生成待对比相似性的新曲线 ,参数1（lineb）进行合并操作
    如果无联合合并线，则返回初始线。
    参数：生成新的非闭合曲线的范围（min_dist）"""
    line_b=line_b0.geometry
    line_a=line_a0.geometry

    gcla1 = 1
    gcla2= 3

    ele = line_b0['ele']
    myid=line_b0['myid']  

    candidates = contour_gdf[((contour_gdf['gcla'] == gcla1)  |(contour_gdf['gcla'] == gcla2)) &  (contour_gdf['ele'] == ele)&(contour_gdf['myid'] != myid)]

    min_dist =500
    line0 = None #为共同进行合并的线
    for _, row in candidates.iterrows():
        current_dist = row.geometry.distance(line_b.centroid)
        if current_dist < min_dist:
            min_dist = current_dist
            line0 = row.geometry      
    if line0 == None:   #未找到联合，进行合并的线,则返回本身 
        return line_b,line_a
    #空间约束
    polygon=line_a.convex_hull
    if not polygon.intersects(line0 ):
        return line_b, line_a

    #依据line0(联合合并线),生成新的闭合线（截取闭合线）
    d = line_b.centroid.distance(line0)  # 中心到联合线的最短距离
    buffer_zone = line0.buffer(d)
    new_line_b = line_b.difference(buffer_zone)
    if new_line_b.geom_type == 'MultiLineString':
        new_line_b = linemerge(new_line_b)
        if new_line_b.geom_type == 'MultiLineString':
            new_line_b = max(new_line_b.geoms, key=lambda x: x.length)

    # 依据新生成的闭合线（lineb），对开放/闭合线（linea）进行截取
    d_buffer1=line_b.centroid.distance(line_a)
    d_buffer2=line_b.centroid.distance(line0)
    d_buffer=max(d_buffer1,d_buffer2)

    buffer_zone = line_b.centroid.buffer(2*d_buffer) #缓冲距离或d
    intersected = line_a.intersection(buffer_zone)

    if isinstance(intersected, MultiLineString):
        new_line_a = max(intersected.geoms, key=lambda x: x.length)
    else:
        new_line_a = intersected

    return new_line_b,new_line_a

# 定义计算混合测度的函数;输入为geose
def cal_fixFea(line):
    # 获取曲线的坐标
    x, y = line.xy
    x = np.array(x)
    y = np.array(y) 
    dx = np.gradient(x)
    dy = np.gradient(y)  

    # 计算角度
    angle = np.arctan2(dy, dx)
    angle_degrees = (np.degrees(angle) + 360) % 360
    fixFea = list(angle_degrees)   
    # # 计算正切值
    # tan_angle = dy / dx
    # fixFea = list(tan_angle)   
    return fixFea

def get_segment_length(line: LineString, min_idx, max_idx) :
    # 依据坐标索引，获取线条长度
    points = list(line.coords)
    sub_points = points[min_idx : max_idx+1]
    
    # 计算相邻点之间的欧氏距离之和
    total_length = 0.0
    for i in range(len(sub_points) -1):
        x1, y1 = sub_points[i]
        x2, y2 = sub_points[i + 1]
        dx, dy = x2 - x1, y2 - y1
        total_length += math.hypot(dx, dy)  # 等同于 sqrt(dx² + dy²)
    
    return total_length


def LCSS_sim(line_an, line_bn):
    """ 连续的公共子序列；tolerance: 最长公共子序列匹配的容忍度    """
    # 计算两个特征序列
    list1 = cal_fixFea(line_an)
    list2 = cal_fixFea(line_bn)
    std_dev1 = np.std(list1)
    std_dev2 = np.std(list2)
    # tolerance=min(std_dev1,std_dev2)
    tolerance=10
    m = len(list1)
    n = len(list2) #闭合线的序列数,即进行了合并操作的线

# DP表和路径追踪矩阵
    dp = [[0]*(n+1) for _ in range(m+1)]
    path = [[[] for _ in range(n+1)] for _ in range(m+1)]  # 记录路径索引对

    for i in range(1, m+1):
        for j in range(1, n+1):
            a, b = list1[i-1], list2[j-1]
            if abs(a - b) < tolerance:
                dp[i][j] = dp[i-1][j-1] + 1
                path[i][j] = path[i-1][j-1] + [(i-1, j-1)]  # 记录原始索引
            else:
                if dp[i-1][j] > dp[i][j-1]:
                    dp[i][j] = dp[i-1][j]
                    path[i][j] = path[i-1][j]
                else:
                    dp[i][j] = dp[i][j-1]
                    path[i][j] = path[i][j-1]

    # 获取最长公共子序列索引对
    lcs_pairs = path[m][n]
    
    # 提取索引范围
    if lcs_pairs:
        list1_indices = [p[0] for p in lcs_pairs]
        list2_indices = [p[1] for p in lcs_pairs]
        min_idx1 = min(list1_indices)
        min_idx2 = min(list2_indices)
        max_idx1 = max(list1_indices)
        max_idx2 = max(list2_indices)
    else:
        min_idx1=max_idx1 = 0
        min_idx2=max_idx2 = 0 # 无匹配时的默认值
    
    lengtha=get_segment_length(line_an,min_idx1,max_idx1)
    lengthb=get_segment_length(line_bn,min_idx2,max_idx2)
    
    # 处理不存在公共序列的情况
    if min(lengtha,lengthb)>0:
        similarity_score=min(lengtha,lengthb)/max(lengtha,lengthb)
    else:
         similarity_score=0

    return similarity_score
################################################################################################################################################################################

""" 闭合-闭合匹配。参数：容差大小和缓冲区距离 """

def merging_CCline(contour_gdf, line_b0, line_a0):
    """ 生成待对比相似性的新曲线 ,lineb需为闭合线;
    如果无联合合并线，则返回初始线;
    参数：联合线的检索范围，高效率（min_dist）"""
    line_b=line_b0.geometry
    line_a=line_a0.geometry

    gcla1= 2
    gcla2= 3

    ele = line_b0['ele']
    myid=line_b0['myid']  

    candidates = contour_gdf[((contour_gdf['gcla'] == gcla1)  |(contour_gdf['gcla'] == gcla2)) &  (contour_gdf['ele'] == ele)&(contour_gdf['myid'] != myid)]

    min_dist = 500 # 减少搜寻时间
    line0 = None #为共同进行合并的线
    for _, row in candidates.iterrows():
        current_dist = row.geometry.distance(line_b.centroid)
        if current_dist < min_dist:
            min_dist = current_dist
            line0 = row.geometry   
            # print(row.myid)   
    if line0 == None:   #未找到联合，进行合并的线,则返回本身 
        return line_b, line_a
    # 增加空间约束
    polygon=line_a.convex_hull
    if not polygon.intersects(line0 ):
        return line_b, line_a

    #依据line0(联合合并线),生成新的闭合线（截取闭合线）
    d = line_b.centroid.distance(line0)  # 中心到联合线的最短距离
    buffer_zone = line0.buffer(d)
    new_line_b = line_b.difference(buffer_zone)
    if new_line_b.geom_type == 'MultiLineString':
        new_line_b = linemerge(new_line_b)
        if new_line_b.geom_type == 'MultiLineString':
            new_line_b = max(new_line_b.geoms, key=lambda x: x.length)

    # 依据新生成的闭合线（lineb），对开放/闭合线（linea）进行截取
    d_buffer1=line_b.centroid.distance(line_a)
    d_buffer2=line_b.centroid.distance(line0)
    d_buffer=max(d_buffer1,d_buffer2)

    buffer_zone = line_b.centroid.buffer(2*d_buffer) #缓冲距离或3d
    intersected = line_a.intersection(buffer_zone)
    
    # # 增加面积限制
    # merged_points = np.vstack([np.array(line_b.coords), np.array(line0.coords)])
    # hullradio=ConvexHull(merged_points).volume/line_a.convex_hull.area
    # if hullradio> 1.33 or hullradio<0.75 :
    #     return line_b, line_a

    if isinstance(intersected, MultiLineString):
        new_line_a = max(intersected.geoms, key=lambda x: x.length)
    else:
        new_line_a = intersected

    return new_line_b,new_line_a


################################################################################################################################################################################
""" 用于截取式的相似性对比 """

def ge_curvesTO(line10, line20):
    """ 获取截断式等高线之间的最相似部分,第一个为截断式,第二个为开放式 """
    line1=line10.geometry
    line2=line20.geometry

    # 检查输入是否为LineString
    if not isinstance(line1, LineString) or not isinstance(line2, LineString):
        raise ValueError("输入必须是LineString类型")

    ##截取line2

    s1 = Point(line1.coords[0])
    e1 = Point(line1.coords[-1])
    p1 = nearest_points(s1, line2)[1]
    p2 = nearest_points(e1, line2)[1]
    

    d1 = line2.project(p1)
    d2 = line2.project(p2)
    start, end = sorted([d1, d2])
    seg2 = substring(line2, start, end)

    ##获取截断的曲线段
    # 在line1上找到最近点p3和p4
    s2 = Point(seg2.coords[0])
    e2 = Point(seg2.coords[-1])   
    p3 = nearest_points(s2, line1)[1]
    p4 = nearest_points(e2, line1)[1]
    
    # 截取line1上p3到p4的曲线段
    d3 = line1.project(p3)
    d4 = line1.project(p4)
    start_seg, end_seg = sorted([d3, d4])
    seg1 = substring(line1, start_seg, end_seg)
    
    seg1_gdf = gpd.GeoDataFrame([line10], geometry=[seg1])
    seg2_gdf = gpd.GeoDataFrame([line20], geometry=[seg2])
    
    return seg1_gdf, seg2_gdf


def ge_curvesTC(line10, line20):
    """ 获取截断式等高线之间的最相似部分,第一个为截断式,第二个为闭合式 """
    line1=line10.geometry
    line2=line20.geometry
    # 检查输入是否为LineString
    if not isinstance(line1, LineString) or not isinstance(line2, LineString):
        raise ValueError("输入必须是LineString类型")

    if line2.is_closed:
        # 在line2上找到最近点p1和p2
        s1, e1 = Point(line1.coords[0]), Point(line1.coords[-1])
        p1 = nearest_points(s1, line2)[1]
        p2 = nearest_points(e1, line2)[1]
        

        # 计算的投影距离
        d1 = line2.project(p1)
        d2 = line2.project(p2)
        start, end = sorted([d1 , d2 ])  

        # 返回两部分
        seg2a = substring(line2, start, end)
        buffer = seg2a.buffer(0.01, quad_segs=8)#构建缓冲区
        seg2b = line2.difference(buffer)
        # print(type(seg2b))
        if isinstance(seg2b, MultiLineString):
            seg2b = linemerge(seg2b)  # 合并相邻线段
            if isinstance(seg2b, MultiLineString):
                seg2b = seg2b.geoms[0]

        ##获取截断的曲线段
        # 在line1上找到最近点p3和p4

        p3 = nearest_points(p1, line1)[1]
        p4 = nearest_points(p2, line1)[1]
        
        # 截取line1上p3到p4的曲线段
        d3 = line1.project(p3)
        d4 = line1.project(p4)
        start_seg, end_seg = sorted([d3, d4])
        seg1 = substring(line1, start_seg, end_seg)
        
        # 计算Hausdorff距离
        hausdorff1 = hausdorff_distance(seg1, seg2a)
        hausdorff2 = hausdorff_distance(seg1, seg2b)
        if hausdorff1 <= hausdorff2:
            seg2=seg2a  
        else:
            seg2=seg2b
    seg1_gdf = gpd.GeoDataFrame([line10], geometry=[seg1])
    seg2_gdf = gpd.GeoDataFrame([line20], geometry=[seg2]) 
    return seg1_gdf, seg2_gdf

def ge_curvesTT(line10, line20):
    """ 获取截断式等高线之间的最相似部分,第1、2个为截断式"""
    line1=line10.geometry
    line2=line20.geometry
    # 检查输入是否为LineString
    if not isinstance(line1, LineString) or not isinstance(line2, LineString):
        raise ValueError("输入必须是LineString类型")

    ##获取截断的曲线段
    # 在line2上找到最近点p1和p2
    s1 = Point(line1.coords[0])
    e1 = Point(line1.coords[-1])
    p1 = nearest_points(s1, line2)[1]
    p2 = nearest_points(e1, line2)[1]
    
    # 截取line2上p1到p2的曲线段
    d1 = line2.project(p1)
    d2 = line2.project(p2)
    start, end = sorted([d1, d2])
    seg2 = substring(line2, start, end)

    ##获取截断的曲线段
    # 在line1上找到最近点p3和p4
    s2 = Point(line2.coords[0])
    e2 = Point(line2.coords[-1])   
    p3 = nearest_points(s2, line1)[1]
    p4 = nearest_points(e2, line1)[1]
    
    # 截取line1上p3到p4的曲线段
    d3 = line1.project(p3)
    d4 = line1.project(p4)
    start_seg, end_seg = sorted([d3, d4])
    seg1 = substring(line1, start_seg, end_seg)
    seg1_gdf = gpd.GeoDataFrame([line10], geometry=[seg1])
    seg2_gdf = gpd.GeoDataFrame([line20], geometry=[seg2])
    return seg1_gdf, seg2_gdf

##################################################################################################

#计算精度

def calculate_metrics4_df(pred, y):
    """ 计算指标值 """
    # 将 DataFrame 列转换为 NumPy 数组
    pred = np.array(pred)
    y = np.array(y)

    # 计算混淆矩阵的各项
    TP = np.sum((pred == 1) & (y == 1))
    TN = np.sum((pred == 0) & (y == 0))
    FP = np.sum((pred == 1) & (y == 0))
    FN = np.sum((pred == 0) & (y == 1))
        
    # 计算准确率、精确率、召回率和 F1 分数
    acc1 = (TP + TN) / (TP + TN + FP + FN)
    precision1 = TP / (TP + FP) if (TP + FP) > 0 else 0
    recall1 = TP / (TP + FN) if (TP + FN) > 0 else 0
    f11 = 2 * precision1 * recall1 / (precision1 + recall1) if (precision1 + recall1) > 0 else 0
    return acc1, precision1, recall1, f11